clc;
clear;

example = 5;
switch example
    case 1
        node = [-1, -1;
            0, -1;
            1, -1;
            -1, 0;
            0, 0;
            1, 0;
            -1, 1;
            0, 1;
            1, 1];
        elem = [5, 4, 1;
            1, 2, 5;
            6, 5, 2;
            2, 3, 6;
            8, 7 ,4;
            4, 5, 8;
            9, 8, 5;
            5, 6, 9];
    case 2
        node = [0, -1;
            -0.707106781187, -0.707106781187;
            0.707106781187, -0.707106781187;
            0,-0.333333334358;
            -1,0;
            -0.707106781187, 0.707106781187;
            1, 0;
            0, 0.333333332309;
            0, 1;
            0.707106781187, 0.707106781187];
        elem = [8, 10,9;
            3, 7, 4;
            4, 8, 5;
            4, 2, 1;
            5, 2, 4;
            4, 1, 3;
            10, 8, 7;
            8, 6, 5;
            7, 8, 4;
            9, 6, 8];
    case 3
        [node, elem] = process_square(0, 1, 0, 1, 3, 3);
    case 4
        r = 1;
        V = [r, 0;
            r/2, (sqrt(3)/2)*r;
            -r/2, (sqrt(3)/2)*r;
            -r, 0;
            -r/2, -(sqrt(3)/2)*r;
            r/2, -(sqrt(3)/2)*r];
        tri = delaunayTriangulation(V);
        node = tri.Points;
        elem = tri.ConnectivityList;
    case 5
        filename = "mesh_hole.msh";
        fileID = fopen(filename, "r");
        data = textscan(fileID, "%f %f %f %f", "Delimiter", " ", "MultipleDelimsAsOne", true);
        fclose(fileID);
        Np = data{1}(1);
        Ne = data{2}(1);
        Nbe = data{3}(1);
        node = [data{1}(2:Np+1), data{2}(2:Np+1)];
        node_out = data{3}(2:Np+1);
        elem = [data{1}(Np+2:Np+Ne+1), data{2}(Np+2:Np+Ne+1), data{3}(Np+2:Np+Ne+1)];
        elem_region = data{4}(Np+2:Np+Ne+1);
        bde_idx = [data{1}(Np+Ne+2:end), data{2}(Np+Ne+2:end), data{3}(Np+Ne+2:end)];
        % case 5 有问题
        % 如果读取网格数据则哟用bde_idx而不是process_refine的bde
        % new_bde生成有问题，导致bdn也有问题
end


basis_type = "P2";
[new_node, new_elem, edge, new_bde, new_bdn] = process_refine(node, elem, basis_type);
mesh.basis_type = basis_type;
mesh.Np = size(node,1);
mesh.Ne = size(elem,1);
mesh.node = new_node;
mesh.Nb = size(new_node,1);
mesh.elem = new_elem;
mesh.edge = edge;
mesh.Ned = size(edge,1);
mesh.bde = new_bde;
mesh.bdn = new_bdn;
process_showmesh(mesh,1);








%% process_square
function [node, elem] = process_square(left, right, bottom, top, Nx, Ny)
[X, Y] = meshgrid(linspace(left, right, Nx+1), linspace(bottom, top, Ny+1));
node = [X(:), Y(:)];
Np = size(node,1);
k = (1:Np-Ny-1)';
cut = (Ny+1)*(1:Nx);
k(cut) = [];
elem = [k, k+Ny+1, k+1;
    k+1, k+Ny+1, k+Ny+2];
elem = sortrows(elem);
end

%% process_refine
function [new_node, new_elem, edge, new_bde, new_bdn] = process_refine(node, elem, basis_type)
Np = size(node,1);
Ne = size(elem,1);
Ned = Np+Ne-1;
elem_ed = zeros(Ne,3);
edge = zeros(Ned,2);
edge_el = edge;
for e = 1:Ne
    n1 = elem(e, 1);
    n2 = elem(e, 2);
    n3 = elem(e, 3);
    tempe = [n1, n2;
        n2, n3;
        n3, n1];
    if e == 1
        edge(1:3,:) = tempe;
        count = 3;
        elem_ed(e,:) = [1, 2, 3];
        edge_el(1:3,1) = 1;
    else
        count1 = 0;
        for m = 1:3
            flag = 1;
            for i = 1:count
                if tempe(m,:) == [edge(i,2), edge(i,1)]
                    flag = 0;
                    elem_ed(e,m) = -i;
                    edge_el(i,2) = e;
                    break;
                end
            end
            if flag == 1
                count1 = count1 + 1;
                new_idx = count + count1;
                edge(new_idx,:) = tempe(m,:);
                elem_ed(e,m) = new_idx;
                edge_el(new_idx,1) = e;
            end
        end
        count = count + count1;
    end
end
bde = [edge_el(edge_el(:,2)==0,1), edge(edge_el(:,2)==0,:), find(edge_el(:,2) == 0)];
keyMap = containers.Map("KeyType", "double", "ValueType", "any");
for i = 1:size(bde,1)
    keyMap(bde(i,2)) = i;
end
new_bde = zeros(size(bde));
new_bde(1,:) = bde(1,:);
for i = 2:size(new_bde,1)
    idx = keyMap(new_bde(i-1,3));
    new_bde(i,:) = bde(idx,:);
end
bdn = unique(new_bde(:,2:3), "stable");
switch basis_type
    case "P1"
        new_node = node;
        new_elem = elem;
        new_bdn = bdn;
    case "P1b"
        center = zeros(Ne,2);
        for e = 1:Ne
            V = [node(elem(e,:),:)];
            center(e,:) = mean(V);
        end
        new_col = (Np+1:Np+Ne)';
        new_node = [node; center];
        new_elem = [elem, new_col];
        new_bdn = bdn;
    case "P2"
        midpoint = 0.5*[(node(edge(:,1),1) + node(edge(:,2),1)), (node(edge(:,1),2) + node(edge(:,2),2))];
        new_node = [node; midpoint];
        new_elem = [elem, abs(elem_ed)+Np];
        new_bdn = reshape([bdn'; (new_bde(:,4)+Np)'], [], 1);
    case "P3"
        differences = node(edge(:,2),:) - node(edge(:,1),:);
        point1_3 = node(edge(:,1),:) + differences/3;
        point2_3 = node(edge(:,1),:) + 2*differences/3;
        triplepoint = [point1_3; point2_3];
        center = zeros(Ne,2);
        for e = 1:Ne
            V = [node(elem(e,:),:)];
            center(e,:) = mean(V);
        end
        new_node = [node; triplepoint; center];
        e1 = elem_ed;
        e1(e1 > 0) = e1(e1 > 0) + Np;
        e1(e1 < 0) = -e1(e1 < 0) + Np+Ned;
        e2 = elem_ed;
        e2(e2 > 0) = e2(e2 > 0) + Np+Ned;
        e2(e2 < 0) = -e2(e2 < 0) + Np;
        new_elem = [elem, e1, e2, (Np+Ned+Ned+1:Np+Ned+Ned+Ne)'];
        new_bdn = reshape([bdn'; (new_bde(:,4)+Np)'; (new_bde(:,4)+Np+Ned)'], [], 1);
    otherwise
        error("Invalid basis type.");
end
end

%% process_plot
function process_plot(mesh, pde, U, basis_type, varargin)
x = mesh.node(:,1);
y = mesh.node(:,2);
if isempty(varargin)
    figure;
    if strcmp(basis_type, "P1")
        x = x(1:mesh.Np);
        y = y(1:mesh.Np);
        elem_re = mesh.elem(:, 1:3);
    elseif strcmp(basis_type, "P1b")
        mod = [1, 2, 4;
            2, 3, 4;
            3, 1, 4];
        elem_re = reshape(mesh.elem(:, mod(:)), 3*mesh.Ne, 3);
    elseif strcmp(basis_type, "P2")
        mod = [1, 4, 6;
            2, 5, 4;
            3, 6, 5;
            4, 5, 6];
        elem_re = reshape(mesh.elem(:, mod(:)), 4*mesh.Ne, 3);
    elseif strcmp(basis_type, "P3")
        mod = [1, 4, 9;
            2, 5, 7;
            3, 6, 8;
            4, 7, 10;
            5, 8, 10;
            6, 9, 10;
            7, 5, 10;
            8, 6, 10;
            9, 4, 10];
        elem_re = reshape(mesh.elem(:, mod(:)), 9*mesh.Ne, 3);
    end
    trisurf(elem_re, x, y, U);
else
    figure;
    scatter3(x, y, U, 20, U, 'filled');
end
colormap('parula');
colorbar;
title(['Example: ', num2str(pde.index)]);
xlabel('x');
ylabel('y');
zlabel('uh');
axis("equal", "tight");
end